setwd("/Users/fernandolucasruiz/Library/CloudStorage/OneDrive-UNIVERSIDADDEMURCIA/Documentos/Fernando/Master Bioinformatica/TFM")
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.2
library(Hmisc)
## Warning: package 'Hmisc' was built under R version 4.3.2
library(limma)
library(AnnotationDbi)
## Warning: package 'AnnotationDbi' was built under R version 4.3.1
## Warning: package 'IRanges' was built under R version 4.3.1
## Warning: package 'S4Vectors' was built under R version 4.3.1
library(org.Hs.eg.db)

library(VennDiagram)
library(ggplot2)
library(gridExtra)
library(purrr)
library(patchwork)
library(nortest)

library(ggrepel)
library(Rtsne)
library(umap)

Leyendo datos

Descargar los datos

# covariables
ROSMAP_RINPMIAGESEX_covs <- readRDS("~/Library/CloudStorage/OneDrive-UNIVERSIDADDEMURCIA/Documentos/Fernando/Master Bioinformatica/TFM/datos/ROSMAP_RINPMIAGESEX_covs.rds")
covs <- ROSMAP_RINPMIAGESEX_covs

#datos corregidos 
ROSMAP_RINPMIAGESEX_resids <- readRDS("~/Library/CloudStorage/OneDrive-UNIVERSIDADDEMURCIA/Documentos/Fernando/Master Bioinformatica/TFM/datos/ROSMAP_RINPMIAGESEX_resids.rds")
mat_exp <- ROSMAP_RINPMIAGESEX_resids
covs$projid <- as.character(covs$projid)
covs$study <- as.factor(covs$study)
covs$ceradsc <- as.factor(covs$ceradsc)
covs$cogdx <- as.factor(covs$cogdx)
covs$neuroStatus <- as.factor(covs$neuroStatus)
rownames(covs) <- covs$mrna_id
describe(covs)
## covs 
## 
##  13  Variables      342  Observations
## --------------------------------------------------------------------------------
## mrna_id 
##        n  missing distinct 
##      342        0      342 
## 
## lowest : 02_120405_2  04_120405_1  05_120405_4  07_120410_4  08_120410_5 
## highest: 950_131107_8 951_131107_8 952_131107_8 954_131107_7 958_131107_7
## --------------------------------------------------------------------------------
## projid 
##        n  missing distinct 
##      342        0      342 
## 
## lowest : 10100574 10101327 10101589 10208143 10222853
## highest: 94722642 94828877 95491648 9841821  98953007
## --------------------------------------------------------------------------------
## study 
##        n  missing distinct 
##      342        0        2 
##                       
## Value        MAP   ROS
## Frequency    167   175
## Proportion 0.488 0.512
## --------------------------------------------------------------------------------
## age_death 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      342        0       20    0.902    86.52    4.746    76.05    79.00 
##      .25      .50      .75      .90      .95 
##    85.00    89.00    90.00    90.00    90.00 
##                                                                             
## Value         67    72    73    74    75    76    77    78    79    80    81
## Frequency      1     6     3     1     2     5    10     5     6     5    10
## Proportion 0.003 0.018 0.009 0.003 0.006 0.015 0.029 0.015 0.018 0.015 0.029
##                                                                 
## Value         82    83    84    85    86    87    88    89    90
## Frequency      8    11    10    20    16    22    20    24   157
## Proportion 0.023 0.032 0.029 0.058 0.047 0.064 0.058 0.070 0.459
## 
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## msex 
##        n  missing distinct 
##      342        0        2 
##                       
## Value          0     1
## Frequency    223   119
## Proportion 0.652 0.348
## --------------------------------------------------------------------------------
## batch 
##        n  missing distinct 
##      342        0        9 
##                                                                 
## Value          0     1     2     3     4     5     6     7     8
## Frequency     11    55    48    47    44    54    38    30    15
## Proportion 0.032 0.161 0.140 0.137 0.129 0.158 0.111 0.088 0.044
## --------------------------------------------------------------------------------
## pmi 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      342        0      196        1    7.378    4.496    3.000    3.557 
##      .25      .50      .75      .90      .95 
##    4.454    6.000    8.562   13.663   18.167 
## 
## lowest : 1       1.33333 1.75    1.83333 1.95   
## highest: 21.4167 22.25   23.3333 23.9167 26     
## --------------------------------------------------------------------------------
## rin 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      342        0       42    0.999    7.049    1.149    5.300    5.600 
##      .25      .50      .75      .90      .95 
##    6.300    7.200    7.875    8.300    8.500 
## 
## lowest : 5   5.1 5.2 5.3 5.4, highest: 8.7 8.8 9   9.1 9.2
## --------------------------------------------------------------------------------
## braaksc 
##        n  missing distinct 
##      342        0        7 
##                                                     
## Value          0     1     2     3     4     5     6
## Frequency      6    34    29    76    94    96     7
## Proportion 0.018 0.099 0.085 0.222 0.275 0.281 0.020
## --------------------------------------------------------------------------------
## ceradsc 
##        n  missing distinct 
##      342        0        4 
##                                   
## Value          1     2     3     4
## Frequency    108    64    38   132
## Proportion 0.316 0.187 0.111 0.386
## --------------------------------------------------------------------------------
## cogdx 
##        n  missing distinct 
##      342        0        5 
##                                         
## Value          1     2     3     4     5
## Frequency    102    62     4   159    15
## Proportion 0.298 0.181 0.012 0.465 0.044
## --------------------------------------------------------------------------------
## apoe_genotype 
##        n  missing distinct 
##      342        0        6 
##                                               
## Value         22    23    24    33    34    44
## Frequency      3    37     7   208    83     4
## Proportion 0.009 0.108 0.020 0.608 0.243 0.012
## --------------------------------------------------------------------------------
## neuroStatus 
##        n  missing distinct 
##      342        0        2 
##                       
## Value          0     1
## Frequency    168   174
## Proportion 0.491 0.509
## --------------------------------------------------------------------------------

geneset de alzheimer

# geneset de Alzheimer extraidos de KEGG

# library(limma)
# library(AnnotationDbi)
# library(org.Hs.eg.db)

tab <- getGeneKEGGLinks(species="hsa")
tab$Symbol <- mapIds(org.Hs.eg.db, tab$GeneID,
                       column="SYMBOL", keytype="ENTREZID")
## 'select()' returned 1:1 mapping between keys and columns
paths <- getKEGGPathwayNames(species="hsa")
paths[paths$PathwayID=="hsa05010", ]
geneset_alz <- tab$Symbol[tab$PathwayID=="hsa05010"]

Geneset en matriz de expresion

Vamos a seleccionar solamente los genes de Alzheimer de KEGG en nuestro dataset de expresión génica.

mat_exp_alz_genes <- mat_exp[, colnames(mat_exp) %in% geneset_alz]

Visualizamos si se han seleccionado todos los genes y cuantos del geneset de KEGG

#library(VennDiagram)
#library(grid)

venn.plot <- venn.diagram(
  x = list(GenesMatriz = colnames(mat_exp), GenesetAlzheimer = geneset_alz),
  category.names = c("Matrix Genes", "Geneset Alzheimer"),
  filename = NULL,
  output = FALSE,  # Asegura que no se exporte a un archivo
  fill = c("#00CED1", "#E9967A"),
  cex = 1,  # Aumenta el tamaño del texto
  fontface = "bold",
  cat.cex = 1,  # Aumenta el tamaño del texto de las categorías
  cat.fontface = "bold",
  cat.default.pos = "text",
  cat.pos = 25, #posicion de las categoricas
  cat.dist = 0.1, #distancia de las categoricas
  rotation.degree = 0, 
  margin = 0.1, # hacerla más pequeña
  lwd=0.5,
  lty = "dashed", # Estilo de línea discontinua
  edge.col = "grey", # Color de los bordes
  main = "Alzheimer genes in exp matrix",
  main.fontface= "bold", 
  main.cex = 2,
  main.pos = c(0.5, 1)
)

grid.newpage()  # Asegura que el lienzo esté limpio
grid.draw(venn.plot)

Sin escalar

PCA datos matriz sin escalar

Hacemos PCA de nuestra nueva matriz con los genes seleccionados de Alzheimer. Los datos vienen libres de sesgos. Tampoco sería necesario el escalado de los datos ya que todos tienen la misma escala y vienen normalizados de antes.

pca_alz_genes <- prcomp(mat_exp_alz_genes)

El 63% de la varinza se explica en la primera componente principal. Eso en una matriz de expresión es raro, ya que suele ser bastante homogéneo.

pca.covs <-  as.data.frame(pca_alz_genes$x[,1:10])
pca.covs <- merge(pca.covs, as.data.frame(covs), by = "row.names")
rownames(pca.covs) <- pca.covs$Row.names
pca.covs$Row.names <- NULL

plots <- list()

# Itera sobre las columnas de covariables y añadir a la lista
for (i in 11:ncol(pca.covs)) {
  # Crea una gráfica para cada variable
  p <- ggplot(pca.covs, aes_string(x = names(pca.covs)[1], y = names(pca.covs)[2], color = names(pca.covs)[i])) +
    geom_point() +
    labs(title = names(pca.covs)[i],
         x = names(pca.covs)[1], y = names(pca.covs)[2], color = names(pca.covs)[i]) +
    modelr::geom_ref_line(h = 0) +
    modelr::geom_ref_line(v = 0) +
    geom_point() +
    xlab("PC 1 ") +
    ylab("PC 2")

  plots[[colnames(pca.covs)[i]]] <- p
}

grid.arrange(plots$study, plots$age_death, plots$msex, plots$batch, plots$pmi, plots$rin, ncol=2)

grid.arrange(plots$braaksc, plots$ceradsc, plots$cogdx, plots$apoe_genotype, plots$neuroStatus, ncol=2)

### Seleccionar muestras con el doble de la desviación típica

Seleccionamos muestras que esten 2 veces la desviacion tipica para la PC1 y PC2

outlierspc1 <- as.data.frame(pca_alz_genes$x[abs(pca_alz_genes$x[,1]) > 2*(pca_alz_genes$sdev[1]),1])

outlierspc2 <- as.data.frame(pca_alz_genes$x[abs(pca_alz_genes$x[,2]) > 2*(pca_alz_genes$sdev[2]),2])

Las ploteamos

Creamos los grupos de muestras seleccionados de los outliers de PC1 y/o PC2

mPC1 <- mat_exp[rownames(mat_exp) %in% rownames(outlierspc1),]
mPC2 <- mat_exp[rownames(mat_exp) %in% rownames(outlierspc2),]
mPC12 <- mat_exp[(rownames(mat_exp) %in% rownames(outlierspc2)) & (rownames(mat_exp) %in% rownames(outlierspc1)),]

not_in_PC12 <- mat_exp[!((rownames(mat_exp) %in% rownames(outlierspc2)) | (rownames(mat_exp) %in% rownames(outlierspc1))),]
# Asegurarse de que todos son vectores del mismo largo para el dataframe final
max_length <- max(length(rownames(mPC1)), length(rownames(mPC2)), length(rownames(mPC12)))

# Normalizar la longitud de los vectores (en caso de que alguno sea más corto)
rownames_mPC1 <- c(rownames(mPC1), rep("", max_length - length(rownames(mPC1))))
rownames_mPC2 <- c(rownames(mPC2), rep("", max_length - length(rownames(mPC2))))
rownames_mPC12 <- c(rownames(mPC12), rep("", max_length - length(rownames(mPC12))))

final_table <- data.frame(
  mPC1 = rownames_mPC1,
  mPC2 = rownames_mPC2,
  mPC1_and_PC2 = rownames_mPC12
)

print(final_table)
##            mPC1         mPC2 mPC1_and_PC2
## 1  500_120515_2 681_120604_1 500_120515_2
## 2  629_120524_2 380_120503_1 629_120524_2
## 3  507_120515_2 500_120515_2 507_120515_2
## 4  605_120523_2 405_120503_2 271_120430_4
## 5  271_120430_4 629_120524_2 367_120502_4
## 6  367_120502_4 507_120515_2             
## 7               657_120529_3             
## 8               316_120501_3             
## 9               584_120522_4             
## 10              271_120430_4             
## 11              367_120502_4             
## 12              337_120501_5             
## 13              232_120425_5             
## 14              234_120425_6             
## 15              790_130530_7
venn.plot <- venn.diagram(
  x = list(mat_exp = rownames(pca_alz_genes$x), mPC1 = rownames(mPC1), mPC2 = rownames(mPC2)),
  category.names = c("Expression matrix", "mPC1", "mPC2" ),
  filename = NULL,
  output = FALSE,  # Asegura que no se exporte a un archivo
  fill = c("#F08080", "#7FFFD4", "#87CEEB"),
  cex = 1,  # Aumenta el tamaño del texto
  fontface = "bold",
  cat.cex = 1,  # Aumenta el tamaño del texto de las categorías
  cat.fontface = "bold",
  cat.default.pos = "text",
  cat.pos = 0, #posicion de las categoricas
  cat.dist = -0.05, #distancia de las categoricas
  rotation.degree = 0, 
  margin = 0.1, # hacerla más pequeña
  lwd=0.5,
  lty = "dashed", # Estilo de línea discontinua
  edge.col = "grey", # Color de los bordes
  main = "Sample set not scaled PCA",
  main.fontface= "bold", 
  main.cex = 2,
  main.pos = c(0.5, 1)
)

grid.newpage()  # Asegura que el lienzo esté limpio
grid.draw(venn.plot)

Covariables en outliers

Metemos una nueva variable en las covariables dependiendo de si son outliers de PC1 o del PC2 o de ninguno.

rownames(covs) <- covs$mrna_id

for (i in rownames(mat_exp)){
  if (i %in% rownames(mPC1) & !(i %in% rownames(mPC12))){
    covs[i, "sampleset_PCA"] <- "mPC1"
  }
  if (i %in% rownames(mPC2) & !(i %in% rownames(mPC12))){
    covs[i, "sampleset_PCA"] <- "mPC2"
  }
  if (i %in% rownames(mPC12)){
    covs[i, "sampleset_PCA"] <- "mPC1 and mPC2"
  }
  if (!(i %in% rownames(mPC1)) & !(i %in% rownames(mPC2))){
    covs[i, "sampleset_PCA"] <- "Not in both"
  }
}

covs$sampleset_PCA <- as.factor(covs$sampleset_PCA)
table(covs$sampleset_PCA)
## 
##          mPC1 mPC1 and mPC2          mPC2   Not in both 
##             1             5            10           326

Ahora vamos a hacer un análisis estadístico para ver si hay diferencias significativas entre las muestras outliers de cada componente principal con las que no lo están (mPC1 vs Not in both, mPC2 vs Not in both, mPC12 vs Not in both)

plots <- list()
# Itera sobre las columnas de covariables y añadir a la lista
for (i in 1:ncol(covs)) {
  col_name <- colnames(covs)[i]
  if (class(covs[[i]]) == "character"){
    next
  }
  if (class(covs[[i]]) == "factor"){
    p <- covs %>%
      group_by(across(all_of(col_name)), sampleset_PCA) %>%
      summarise(count = n(), .groups = 'drop') %>%
      mutate(total = sum(count), percentage = (count / total) * 100) %>%
      ggplot(aes_string(x = names(covs)[i], y = "percentage", fill = "sampleset_PCA")) +
      geom_bar(stat = "identity", position = "fill") +
      labs(x = NULL, y = "%", fill = col_name) +
      theme_minimal()
    
    plots[[colnames(covs)[i]]] <- p
  }
  if (class(covs[[i]]) == "numeric"){
    p <- ggplot(covs, aes_string(y=names(covs)[i], x = "sampleset_PCA", fill = "sampleset_PCA")) +
      geom_boxplot(alpha=.7)+
      geom_jitter(color ="#8B1A1A", size = 0.7, alpha = 0.7) +
      theme_minimal() +
      labs(title = names(covs)[i], x = names(covs)[i])
    
    plots[[colnames(covs)[i]]] <- p
  }
}

grid.arrange(plots$age_death, plots$pmi, plots$rin, ncol=2)

grid.arrange(plots$study , plots$msex, plots$batch, plots$braaksc, ncol=2)

grid.arrange(plots$ceradsc, plots$cogdx, plots$apoe_genotype,  plots$neuroStatus, ncol=2)

statistical.tests <- function(df, outlier){
  for (i in colnames(df)) {
    #outlier <- "mPC1" 
      if (i == "mrna_id" | i == "projid" | i == "sampleset_PCA" | i == "sampleset_tSNE")  {
          next
      }
    
    # filtro primero los datos por los que quiero comparar outlier y no incluidos
      filtered_df <- df[df$sampleset_PCA %in% c(outlier, "Not in both"), ]
    
      # si la variable es numérica 
      if (class(filtered_df[[i]]) == "numeric") {
        
        # hacemos test de normalidad de la variable para ver si hacemos test paramétrico o no paramétrico
        normality <- lillie.test(filtered_df[[i]])
        
        # si los datos no son paramétricos
        if (normality$p.value < 0.05) { 
          
          # Verificar que sampleset_PCA tenga exactamente 2 niveles después del filtrado
          if (length(unique(filtered_df$sampleset_PCA)) == 2) {
              wilcox_test <- wilcox.test(filtered_df[[i]] ~ filtered_df$sampleset_PCA)
              
              if (wilcox_test$p.value < 0.05) {
                  print(paste("U Mann-Witney test significativo de", i, "con un pvalor de:", round(wilcox_test$p.value, 6), "entre", outlier, " y 'Not in both'"))
                
              } else {
                print(paste("No hay diferencia significativa de la variable", i, "entre", outlier, "y 'Not in both'"))
              }
          }
          
          # si la variable numérica es paramétrica
         } else if (normality$p.value >= 0.05) {
  
            # hacemos test de la varianza
            var_test <- var.test(filtered_df[[i]] ~ filtered_df$sampleset_PCA)
          
            # le pasamos el t.test dependiendo de si las varanzas son iguales o no
            if (var_test$p.value < 0.05) {
              t_test <-  t.test(filtered_df[[i]] ~ filtered_df$sampleset_PCA, var.equal = F)
              if (t_test$p.value < 0.05){
                print(paste("T test significativo (varianzas desiguales) de",
                            i, "con un pvalor de:", round(t_test$p.value, 6)))
              } else {
                print(paste("No hay diferencia significativa de la variable", i,  "entre", outlier, "y 'Not in both'"))
              }
              
            } else {
              t_test <-  t.test(filtered_df[[i]] ~ filtered_df$sampleset_PCA, var.equal = T)
            
              if (t_test$p.value < 0.05){
                print(paste("T test significativo (varianzas iguales) de",
                          i, "con un pvalor de:", round(t_test$p.value, 6)))
              } else {
                print(paste("No hay diferencia significativa de la variable", i,  "entre", outlier, "y 'Not in both'"))
              }
            }
         } 
      } 
    
      # si la variable es categórica le vamos a pasar un test chisq
      if (class(filtered_df[[i]]) == "factor") {
        
        # hacemos la tabla de contingencia
        contingency_table <- table(filtered_df[[i]], filtered_df$sampleset_PCA)
        
        # test de Chi-cuadrado
        test_chi <- chisq.test(contingency_table)
        
        if (!is.na(test_chi$p.value) && test_chi$p.value < 0.05) {
            print(paste("Test de Chi-cuadrado significativo para", i, "con un p-valor de:", round(test_chi$p.value, 6), "entre", outlier, " y 'Not in both'"))
        } else {
            print(paste("No hay diferencia significativa de la variable", i,  "entre", outlier, "y 'Not in both'"))
          }
      }
  }
}
statistical.tests(covs, "mPC1")
## [1] "No hay diferencia significativa de la variable study entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable age_death entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable pmi entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable rin entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC1 y 'Not in both'"
statistical.tests(covs, "mPC2")
## [1] "No hay diferencia significativa de la variable study entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable age_death entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable pmi entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable rin entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC2 y 'Not in both'"
statistical.tests(covs, "mPC12")
## [1] "No hay diferencia significativa de la variable study entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC12 y 'Not in both'"

Escalado

PCA datos matriz escalando

pca_alz_genes_scaled <- prcomp(mat_exp_alz_genes, scale = T)
# Scree plot con los datos escalados
var_exp <- pca_alz_genes_scaled$sdev^2
prop_var_exp <- var_exp / sum(var_exp)
cum_var_exp <- cumsum(prop_var_exp)

df_var_exp <- data.frame(Comp = 1:length(prop_var_exp), VarExp = prop_var_exp)
df_cum_var_exp <- data.frame(Comp = 1:length(cum_var_exp), CumVarExp = cum_var_exp)

ggplot(df_var_exp[1:20,], aes(x = Comp, y = VarExp)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_line(aes(group = 1), color = "blue") +
  geom_point(color = "blue") +
  theme_minimal() +
  labs(x = "Principal components", y = "Variance", title = "Scree Plot scaled") +
  ylim(c(0,1)) +
  geom_line(data = df_cum_var_exp[1:20,], aes(x = Comp, y = CumVarExp), color="#8B1A1A") +
  geom_point(data = df_cum_var_exp[1:20,], aes(x = Comp, y = CumVarExp), color = "red") +
  geom_bar(data = df_cum_var_exp[1:20,], aes(x = Comp, y = CumVarExp), stat = "identity", fill = "red", alpha= 0.25) +
  annotate("text", x = 4, y = 0.85, label = "Cumulative Scree Plot", color = "#8B1A1A", size = 4) +
  geom_text(data = df_cum_var_exp[1:20,], aes(x=Comp, y = CumVarExp +0.04, label = round(CumVarExp, 2)))

Vemos que con las 20 primeras no encontramos ni siquiera el 80% de la varianza.

pca.scaled.covs <-  as.data.frame(pca_alz_genes_scaled$x[,1:10])
pca.scaled.covs <- merge(pca.scaled.covs, as.data.frame(covs), by = "row.names")
rownames(pca.scaled.covs) <- pca.scaled.covs$Row.names
pca.scaled.covs$Row.names <- NULL

plots <- list()

# Itera sobre las columnas de covariables y añadir a la lista
for (i in 11:ncol(pca.scaled.covs)) {
  # Crea una gráfica para cada variable
  p <- ggplot(pca.scaled.covs, aes_string(x = names(pca.scaled.covs)[1], y = names(pca.scaled.covs)[2], color = names(pca.scaled.covs)[i])) +
    geom_point() +
    labs(title = names(pca.scaled.covs)[i],
         x = names(pca.scaled.covs)[1], y = names(pca.scaled.covs)[2], color = names(pca.scaled.covs)[i]) +
    modelr::geom_ref_line(h = 0) +
    modelr::geom_ref_line(v = 0) +
    geom_point() +
    xlab("PC 1 ") +
    ylab("PC 2")

  plots[[colnames(pca.scaled.covs)[i]]] <- p
}

grid.arrange(plots$study, plots$age_death, plots$msex, plots$batch, plots$pmi, plots$rin, ncol=2)

grid.arrange(plots$braaksc, plots$ceradsc, plots$cogdx, plots$apoe_genotype, plots$neuroStatus, ncol=2)

Seleccionar muestras con el doble de la desviación típica

outlierspc1_scaled <- as.data.frame(pca_alz_genes_scaled$x[abs(pca_alz_genes_scaled$x[,1]) > 2*(pca_alz_genes_scaled$sdev[1]),1])

outlierspc2_scaled <- as.data.frame(pca_alz_genes_scaled$x[abs(pca_alz_genes_scaled$x[,2]) > 2*(pca_alz_genes_scaled$sdev[2]),2])
df <- as.data.frame(pca_alz_genes_scaled$x)

plot1 <- ggplot(df, aes(x = PC1)) +
  geom_density(fill = "#00CED1", alpha = 0.5) +
  geom_vline(xintercept = mean(df$PC1) + 2*pca_alz_genes_scaled$sdev[1], linetype = "dashed", color = "blue") +
  geom_vline(xintercept = mean(df$PC1) - 2*pca_alz_genes_scaled$sdev[1], linetype = "dashed", color = "blue") +
  labs(title = "PC1 density with 2*SD scaled PCA") +
  geom_text(data = outlierspc1_scaled, 
            aes(x = outlierspc1_scaled[,1], y= 0, label = rownames(outlierspc1_scaled)), 
            vjust = 1.5, hjust = 0, size = 3, color = "#E9965A", angle = 90, fontface= "bold") +
  geom_rug(data = as.data.frame(pca_alz_genes_scaled$x), aes(x= PC1, y = 0), color= ifelse(abs(pca_alz_genes_scaled$x[,1]) > 2*(pca_alz_genes_scaled$sdev[1]), "#8B1A1A", "grey")) +
  theme_minimal()

plot2 <- ggplot(df, aes(x = PC2)) +
  geom_density(fill = "#00CED1", alpha = 0.5) +
  geom_vline(xintercept = mean(df$PC2) + 2*pca_alz_genes_scaled$sdev[2], linetype = "dashed", color = "blue") +
  geom_vline(xintercept = mean(df$PC2) - 2*pca_alz_genes_scaled$sdev[2], linetype = "dashed", color = "blue") +
  labs(title = "PC2 density with 2*SD scaled PCA") +
  geom_text(data = outlierspc2_scaled, 
            aes(x = outlierspc2_scaled[,1], y= 0, label = rownames(outlierspc2_scaled)), 
            vjust = 1.5, hjust = 0, size = 3, color = "#E9965A", angle = 90, fontface= "bold") +
  geom_rug(data = as.data.frame(pca_alz_genes_scaled$x), aes(x= PC2, y = 0), color= ifelse(abs(pca_alz_genes_scaled$x[,2]) > 2*(pca_alz_genes_scaled$sdev[2]), "#8B1A1A", "grey")) +
  theme_minimal()

grid.arrange(plot1, plot2, ncol = 1)

mPC1_scaled <- mat_exp[rownames(mat_exp) %in% rownames(outlierspc1_scaled),]
mPC2_scaled <- mat_exp[rownames(mat_exp) %in% rownames(outlierspc2_scaled),]
mPC12_scaled <- mat_exp[(rownames(mat_exp) %in% rownames(outlierspc2_scaled)) & (rownames(mat_exp) %in% rownames(outlierspc1_scaled)),]

not_in_PC12_scaled <- mat_exp[!((rownames(mat_exp) %in% rownames(outlierspc2_scaled)) | (rownames(mat_exp) %in% rownames(outlierspc1_scaled))),]
# Asegurarse de que todos son vectores del mismo largo para el dataframe final
max_length <- max(length(rownames(mPC1_scaled)), length(rownames(mPC2_scaled)), length(rownames(mPC12_scaled)))

# Normalizar la longitud de los vectores (en caso de que alguno sea más corto)
rownames_mPC1_scaled <- c(rownames(mPC1_scaled), rep("", max_length - length(rownames(mPC1_scaled))))
rownames_mPC2_scaled <- c(rownames(mPC2_scaled), rep("", max_length - length(rownames(mPC2_scaled))))
rownames_mPC12_scaled <- c(rownames(mPC12_scaled), rep("", max_length - length(rownames(mPC12_scaled))))

final_table <- data.frame(
  mPC1 = rownames_mPC1_scaled,
  mPC2 = rownames_mPC2_scaled,
  mPC1_and_PC2 = rownames_mPC12_scaled
)

print(final_table)
##            mPC1         mPC2 mPC1_and_PC2
## 1  274_120430_1 384_120503_1             
## 2  724_120605_2 380_120503_1             
## 3  534_120516_2 300_120430_2             
## 4  657_120529_3  07_120410_4             
## 5  193_120424_3 186_120424_4             
## 6  584_120522_4 628_120524_4             
## 7  271_120430_4 336_120501_5             
## 8  337_120501_5 691_120605_5             
## 9   79_120417_5  08_120410_5             
## 10 215_120425_5 119_120418_6             
## 11 232_120425_5 893_130923_7             
## 12 505_120515_6 905_131010_7             
## 13 537_120516_6                          
## 14 769_130523_7                          
## 15 929_131031_8
venn.plot <- venn.diagram(
  x = list(mat_exp = rownames(mat_exp), mPC1 = rownames(mPC1_scaled), mPC2 = rownames(mPC2_scaled)),
  category.names = c("Expression matrix", "mPC1", "mPC2" ),
  filename = NULL,
  output = FALSE,  # Asegura que no se exporte a un archivo
  fill = c("#F08080", "#7FFFD4", "#87CEEB"),
  cex = 1,  # Aumenta el tamaño del texto
  fontface = "bold",
  cat.cex = 1,  # Aumenta el tamaño del texto de las categorías
  cat.fontface = "bold",
  cat.default.pos = "text",
  cat.pos = 0, #posicion de las categoricas
  cat.dist = -0.05, #distancia de las categoricas
  rotation.degree = 0, 
  margin = 0.1, # hacerla más pequeña
  lwd=0.5,
  lty = "dashed", # Estilo de línea discontinua
  edge.col = "grey", # Color de los bordes
  main = "Sample set scaled PCA",
  main.fontface= "bold", 
  main.cex = 2,
  main.pos = c(0.5, 1)
)

grid.newpage()  # Asegura que el lienzo esté limpio
grid.draw(venn.plot)

Voy ahora a ver de esas muestras seleccionadas, qué covariables están asociadas con esas muestras outliers

Covariables en outliers

rownames(covs) <- covs$mrna_id
covs.scaled <- covs

for (i in rownames(covs.scaled)){
  if (i %in% rownames(mPC1_scaled) & !(i %in% rownames(mPC12_scaled))){
    covs.scaled[i, "sampleset_PCA"] <- "mPC1"
  }
  if (i %in% rownames(mPC2_scaled) & !(i %in% rownames(mPC12_scaled))){
    covs.scaled[i, "sampleset_PCA"] <- "mPC2"
  }
  if (i %in% rownames(mPC12_scaled)){
    covs.scaled[i, "sampleset_PCA"] <- "mPC1 and mPC2"
  }
  if (!(i %in% rownames(mPC1_scaled)) & !(i %in% rownames(mPC2_scaled))){
    covs.scaled[i, "sampleset_PCA"] <- "Not in both"
  }
}

covs.scaled$sampleset_PCA <- as.factor(covs.scaled$sampleset_PCA)
table(covs.scaled$sampleset_PCA)
## 
##          mPC1 mPC1 and mPC2          mPC2   Not in both 
##            15             0            12           315
plots <- list()
# Itera sobre las columnas de covariables y añadir a la lista
for (i in 1:ncol(covs.scaled)) {
  col_name <- colnames(covs.scaled)[i]
  if (class(covs.scaled[[i]]) == "character"){
    next
  }
  if (class(covs.scaled[[i]]) == "factor"){
    p <- covs.scaled %>%
      group_by(across(all_of(col_name)), sampleset_PCA) %>%
      summarise(count = n(), .groups = 'drop') %>%
      mutate(total = sum(count), percentage = (count / total) * 100) %>%
      ggplot(aes_string(x = names(covs.scaled)[i], y = "percentage", fill = "sampleset_PCA")) +
      geom_bar(stat = "identity", position = "fill") +
      labs(x = NULL, y = "%", fill = col_name) +
      theme_minimal()
    
    plots[[colnames(covs.scaled)[i]]] <- p
  }
  if (class(covs.scaled[[i]]) == "numeric"){
    p <- ggplot(covs.scaled, aes_string(y=names(covs.scaled)[i], x = "sampleset_PCA", fill = "sampleset_PCA")) +
      geom_boxplot(alpha=.7)+
      geom_jitter(color ="#8B1A1A", size = 0.7, alpha = 0.7) +
      theme_minimal() +
      labs(title = names(covs.scaled)[i], x = names(covs.scaled)[i])
    
    plots[[colnames(covs.scaled)[i]]] <- p
  }
}

grid.arrange(plots$age_death, plots$pmi, plots$rin, ncol=2)

grid.arrange(plots$study , plots$msex, plots$batch, plots$braaksc, ncol=2)

grid.arrange(plots$ceradsc, plots$cogdx, plots$apoe_genotype,  plots$neuroStatus, ncol=2)

statistical.tests(covs.scaled, "mPC1")
## [1] "No hay diferencia significativa de la variable study entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable age_death entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable pmi entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable rin entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC1 y 'Not in both'"
statistical.tests(covs.scaled, "mPC2")
## [1] "No hay diferencia significativa de la variable study entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable age_death entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable pmi entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable rin entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC2 y 'Not in both'"
statistical.tests(covs.scaled, "mPC12")
## [1] "No hay diferencia significativa de la variable study entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC12 y 'Not in both'"

tSNE

Parámetros clave en t-SNE: Perplexity (Perplejidad): Este es uno de los parámetros más importantes en t-SNE. La perplejidad se relaciona con el número de vecinos cercanos que t-SNE considera cuando mapea los datos. Un valor demasiado bajo hace que el modelo se concentre demasiado en la estructura local, mientras que un valor demasiado alto puede llevar a una visión global que pierde detalles. La elección óptima depende del tamaño del dataset, pero valores típicos están entre 5 y 50.

Número de iteraciones: t-SNE es un algoritmo iterativo, y el número de iteraciones puede afectar la estabilidad de los resultados. Un número insuficiente de iteraciones puede resultar en una organización incompleta, mientras que demasiadas iteraciones pueden no proporcionar beneficios adicionales y aumentar el tiempo de cálculo.

Tasa de aprendizaje: Este parámetro controla el tamaño del paso en cada actualización durante la optimización. Una tasa de aprendizaje muy baja puede hacer que el algoritmo tarde mucho en converger, mientras que una tasa demasiado alta puede llevar a una convergencia pobre.

set.seed(1234)
tsne <- Rtsne(mat_exp_alz_genes, dims = 2, theta = 0.0, )

tsne.data <- as.data.frame(tsne$Y)
row.names(tsne.data) <- row.names(mat_exp_alz_genes)
tsne.data.covs <- merge(tsne.data, covs, by = "row.names")
tsne.data.covs$Row.names <- NULL
row.names(tsne.data.covs) <- tsne.data.covs$mrna_id

ggplot(tsne.data.covs, aes(x=V1, y=V2, color=neuroStatus)) +  
  geom_point(size=2) +
  guides(colour=guide_legend(override.aes=list(size=6))) +
  xlab("tSNE1") + ylab("tSNE2") +
  ggtitle("t-SNE") +
  theme_minimal() +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank()) +
  scale_colour_brewer(palette = "Set1")

Seleccionamos outliers

rownames(tsne.data) <- rownames(mat_exp_alz_genes)

# Calcular la media y la desviación estándar para cada componente de t-SNE
mean.tsne1 <- mean(tsne.data[,1])
sd.tsne1 <- sd(tsne.data[,1])
mean.tsne2 <- mean(tsne.data[,2])
sd.tsne2 <- sd(tsne.data[,2])

# Identificar muestras a más de 2 desviaciones estándar de la media
outliers.tsne1 <- tsne.data[abs(tsne.data[,1] - mean.tsne1) > 2 * sd.tsne1, ]
outliers.tsne2 <-tsne.data[abs(tsne.data[,2] - mean.tsne2) > 2 * sd.tsne2, ]
## Warning in mean.default(df$V2): argument is not numeric or logical: returning
## NA

## Warning in mean.default(df$V2): argument is not numeric or logical: returning
## NA
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_vline()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_vline()`).

mtSNE1 <- mat_exp[rownames(mat_exp) %in% rownames(outliers.tsne1),]
mtSNE2 <- mat_exp[rownames(mat_exp) %in% rownames(outliers.tsne2),]
mtSNE12 <- mat_exp[(rownames(mat_exp) %in% rownames(outliers.tsne1)) & (rownames(mat_exp) %in% rownames(outliers.tsne2)),]

not_in_tSNE12 <- mat_exp[!((rownames(mat_exp) %in% rownames(outliers.tsne1)) | (rownames(mat_exp) %in% rownames(outliers.tsne2))),]
# Asegurarse de que todos son vectores del mismo largo para el dataframe final
max_length <- max(length(rownames(mtSNE1)), length(rownames(mtSNE2)), length(rownames(mtSNE12)))

# Normalizar la longitud de los vectores (en caso de que alguno sea más corto)
rownames_tsne1 <- c(rownames(mtSNE1), rep("", max_length - length(rownames(mtSNE1))))
rownames_tsne2 <- c(rownames(mtSNE2), rep("", max_length - length(rownames(mtSNE2))))
rownames_tsne12 <- c(rownames(mtSNE12), rep("", max_length - length(rownames(mtSNE12))))

final_table <- data.frame(
  mtSNE1 = rownames_tsne1,
  mtSNE2 = rownames_tsne2,
  mtSNE1_and_mtSNE2 = rownames_tsne12
)

print(final_table)
##          mtSNE1       mtSNE2 mtSNE1_and_mtSNE2
## 1  681_120604_1 666_120530_1                  
## 2  790_130530_7 426_120507_1                  
## 3               430_120507_1                  
## 4               380_120503_1                  
## 5               500_120515_2                  
## 6               629_120524_2                  
## 7               507_120515_2                  
## 8               544_120516_2                  
## 9               203_120424_2                  
## 10              605_120523_2                  
## 11              568_120521_3                  
## 12              292_120430_3                  
## 13              186_120424_4                  
## 14              367_120502_4                  
## 15              327_120501_5                  
## 16              336_120501_5                  
## 17              257_120426_5
venn.plot <- venn.diagram(
  x = list(mat_exp = rownames(tsne.data), mtSNE1 = rownames(mtSNE1), mtSNE2 = rownames(mtSNE2)),
  category.names = c("Expression matrix", "mtSNE1", "mtSNE2" ),
  filename = NULL,
  output = FALSE,  # Asegura que no se exporte a un archivo
  fill = c("#F08080", "#7FFFD4", "#87CEEB"),
  cex = 1,  # Aumenta el tamaño del texto
  fontface = "bold",
  cat.cex = 1,  # Aumenta el tamaño del texto de las categorías
  cat.fontface = "bold",
  cat.default.pos = "text",
  cat.pos = 0, #posicion de las categoricas
  cat.dist = -0.05, #distancia de las categoricas
  rotation.degree = 0, 
  margin = 0.1, # hacerla más pequeña
  lwd=0.5,
  lty = "dashed", # Estilo de línea discontinua
  edge.col = "grey", # Color de los bordes
  main = "Sample set tSNE",
  main.fontface= "bold", 
  main.cex = 2,
  main.pos = c(0.5, 1)
)

grid.newpage()  # Asegura que el lienzo esté limpio
grid.draw(venn.plot)

for (i in rownames(covs)){
  if (i %in% rownames(mtSNE1) & !(i %in% rownames(mtSNE12))){
    covs[i, "sampleset_tSNE"] <- "mtSNE1"
  }
  if (i %in% rownames(mtSNE2) & !(i %in% rownames(mtSNE12))){
    covs[i, "sampleset_tSNE"] <- "mtSNE2"
  }
  if (i %in% rownames(mtSNE12)){
    covs[i, "sampleset_tSNE"] <- "mtSNE1 and mtSNE2"
  }
  if (!(i %in% rownames(mtSNE1)) & !(i %in% rownames(mtSNE2))){
    covs[i, "sampleset_tSNE"] <- "Not in both"
  }
}

covs$sampleset_tSNE <- as.factor(covs$sampleset_tSNE)
table(covs$sampleset_tSNE)
## 
##      mtSNE1      mtSNE2 Not in both 
##           2          17         323
plots <- list()
# Itera sobre las columnas de covariables y añadir a la lista
for (i in 1:ncol(covs)) {
  col_name <- colnames(covs)[i]
  
  if (class(covs[[i]]) == "character"){
    next
  }
  if (class(covs[[i]]) == "factor"){
    p <- covs %>%
      group_by(across(all_of(col_name)), sampleset_tSNE) %>%
      summarise(count = n(), .groups = 'drop') %>%
      mutate(total = sum(count), percentage = (count / total) * 100) %>%
      ggplot(aes_string(x = names(covs)[i], y = "percentage", fill = "sampleset_tSNE")) +
      geom_bar(stat = "identity", position = "fill") +
      labs(x = NULL, y = "%", fill = col_name) +
      theme_minimal()
    
    plots[[colnames(covs)[i]]] <- p
  }
  if (class(covs[[i]]) == "numeric"){
    p <- ggplot(covs, aes_string(y=names(covs)[i], x = "sampleset_tSNE", fill = "sampleset_tSNE")) +
      geom_boxplot(alpha=.7)+
      geom_jitter(color ="#8B1A1A", size = 0.7, alpha = 0.7) +
      theme_minimal() +
      labs(title = names(covs)[i], x = names(covs)[i])
    
    plots[[colnames(covs)[i]]] <- p
  }
}

grid.arrange(plots$age_death, plots$pmi, plots$rin, ncol=2)

grid.arrange(plots$study , plots$msex, plots$batch, plots$braaksc, ncol=2)

grid.arrange(plots$ceradsc, plots$cogdx, plots$apoe_genotype,  plots$neuroStatus, ncol=2)

statistical.tests(covs, "mtSNE1")
## [1] "No hay diferencia significativa de la variable study entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mtSNE1 y 'Not in both'"
statistical.tests(covs, "mtSNE2")
## [1] "No hay diferencia significativa de la variable study entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mtSNE2 y 'Not in both'"
statistical.tests(covs, "mtSNE12")
## [1] "No hay diferencia significativa de la variable study entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mtSNE12 y 'Not in both'"

UMAP

local.config <- umap.defaults
# local.config$n_neighbors <- 4
# local.config$n_components <- 2
# local.config$n_epochs <- 100
# local.config$metric<- "euclidean"
set.seed(1234)
umap.ad <- umap(mat_exp_alz_genes,random_stage=1234, local.config)

umap.data <- as.data.frame(umap.ad$layout)
row.names(umap.data) <- row.names(mat_exp_alz_genes)
umap.data.covs <- merge(umap.data, covs, by = "row.names")
umap.data.covs$Row.names <- NULL
row.names(umap.data.covs) <- umap.data.covs$mrna_id

ggplot(umap.data.covs, aes(x=V1, y=V2, color = neuroStatus)) +  
  geom_point(size=2) +
  guides(colour=guide_legend(override.aes=list(size=6))) +
  xlab("UMAP1") + ylab("UMAP2") +
  ggtitle("UMAP") +
  theme_minimal() +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank()) +
  scale_colour_brewer(palette = "Set2")